%
% Jean-Sébastien Parent & Maxime Dion, summer 2008
%
%

\defmodule {VarianceGammaProcessDiffPCA}

Same as \class{VarianceGammaProcessDiff}, but the two inner 
\class{GammaProcess}'es are of PCA type.  Also, 
\texttt{generatePath(double[] uniforms01)} distributes 
the uniform random variates to the \class{GammaProcessPCA}'s  according to their 
eigenvalues, i.e. the \class{GammaProcessPCA} with the higher eigenvalue 
gets the next uniform random number.  If one should decide to create a 
\class{VarianceGammaProcessDiffPCA} by giving two \class{GammaProcessPCA}'s to an 
objet of the class \class{VarianceGammaProcessDiff}, the uniform random 
numbers would not be given this way
to the \class{GammaProcessPCA}'s; this  might give less variance reduction when 
used with QMC.

\bigskip\hrule\bigskip

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{code}
\begin{hide}
/*
 * Class:        VarianceGammaProcessDiffPCA
 * Description:  
 * Environment:  Java
 * Software:     SSJ 
 * Copyright (C) 2001  Pierre L'Ecuyer and Universite de Montreal
 * Organization: DIRO, Universite de Montreal
 * @author       
 * @since        2008

 * SSJ is free software: you can redistribute it and/or modify it under
 * the terms of the GNU General Public License (GPL) as published by the
 * Free Software Foundation, either version 3 of the License, or
 * any later version.

 * SSJ is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.

 * A copy of the GNU General Public License is available at
   <a href="http://www.gnu.org/licenses">GPL licence site</a>.
 */
\end{hide}
package umontreal.iro.lecuyer.stochprocess;\begin{hide}
import umontreal.iro.lecuyer.rng.*;
import umontreal.iro.lecuyer.probdist.*;
import umontreal.iro.lecuyer.randvar.*;

\end{hide}

public class VarianceGammaProcessDiffPCA extends VarianceGammaProcessDiff \begin{hide} {
    int[] indexEigenUp;
    int[] indexEigenDw;

\end{hide}
\end{code}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Constructors}
\begin{code}

   public VarianceGammaProcessDiffPCA (double s0, double theta,
                                       double sigma, double nu,
                                       RandomStream stream) \begin{hide} {
     super(s0, theta, sigma, nu, 
	  new GammaProcessPCA (0.0, 1.0, 1.0, stream),
	  new GammaProcessPCA (0.0, 1.0, 1.0, stream));
    // Params mu, nu of the 2 gamma processes are redefined in init()
    // which will be called after a call to 'setObservTimes'
}\end{hide}
\end{code}
\begin{tabb} Constructs a new \class{VarianceGammaProcessDiffPCA} with 
parameters  $\theta = \texttt{theta}$, $\sigma = \texttt{sigma}$, $\nu = \texttt{nu}$ 
and initial value $S(t_{0}) = \texttt{s0}$.  There is only
one \externalclass{umontreal.iro.lecuyer.rng}{RandomStream} here which is
used for the two inner \class{GammaProcessPCA}'s.  The other
parameters are set as in \class{VarianceGammaProcessDiff}.
\end{tabb}
\begin{code}

   public VarianceGammaProcessDiffPCA (double s0, double theta,
                                       double sigma, double nu,
                                       GammaProcessPCA gpos,
                                       GammaProcessPCA gneg) \begin{hide} {
    super(s0, theta, sigma, nu, gpos, gneg); // from VarianceGammaProcessDiff
    // Params mu, nu of the 2 gamma processes are redefined in init()
    // which will be called after a call to 'setObservTimes'
}\end{hide}
\end{code}
\begin{tabb} Constructs a new \class{VarianceGammaProcessDiffPCA} with 
parameters  $\theta = \texttt{theta}$, $\sigma = \texttt{sigma}$, $\nu = \texttt{nu}$ 
and initial value $S(t_{0}) = \texttt{s0}$.  As in 
\class{VarianceGammaProcessDiff}, the 
\externalclass{umontreal.iro.lecuyer.rng}{RandomStream} of \texttt{gneg} is replaced by
the one of \texttt{gpos} to avoid any confusion.
\end{tabb}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection* {Methods}
\begin{code} 

   public double nextObservation() \begin{hide} {
        throw new UnsupportedOperationException 
        ("Impossible with PCA, use generatePath() instead.");
    }\end{hide}
\end{code}
\begin{tabb} This method is not implemented is this class since
the path cannot be generated sequentially.
\end{tabb}
\begin{code}
\begin{hide} 

   public double[] generatePath() {
        double[] u = new double[2*d];
        for(int i =0; i < 2*d; i++)
            u[i] = getStream().nextDouble();
        return generatePath(u);
    }

   public double[] generatePath(double[] uniform01)  {
        int dd = uniform01.length;
        int d = dd / 2;

        if(dd % 2 != 0){
            throw new IllegalArgumentException (
                     "The Array uniform01 must have a even length");
        }

        double[] QMCpointsUP = new double[d];
        double[] QMCpointsDW = new double[d];

        for(int i = 0; i < d; i++){
             QMCpointsUP[i] = uniform01[ indexEigenUp[i] ];
             QMCpointsDW[i] = uniform01[ indexEigenDw[i] ];
        }
        gpos.resetStartProcess();
        gneg.resetStartProcess();

        double[] pathUP = gpos.generatePath(QMCpointsUP);
        double[] pathDOWN = gneg.generatePath(QMCpointsDW);

        for (int i=0; i<d; i++) {
           path[i+1] = x0 + pathUP[i+1] - pathDOWN[i+1];
        }
        observationIndex   = d;
        observationCounter = d;
        return path;
    }


   protected void init() {
        super.init ();  // from VarianceGammaProcessDiff
	if( observationTimesSet){
        // Two lines below (casts) should be reinstated after fix inheritance PCA/PCABridge.
	    double[] eigenValUp = ((GammaProcessPCA)gpos).getBMPCA().getSortedEigenvalues();
	    double[] eigenValDw = ((GammaProcessPCA)gneg).getBMPCA().getSortedEigenvalues();
	    indexEigenUp = new int[d];
	    indexEigenDw = new int[d];

	    int iUp = 0;
	    int iDw = 0;
	    for(int iQMC = 0; iQMC < 2*d; iQMC++){
		if(iUp == d) {indexEigenDw[iDw] = iQMC; iDw++;continue;}
        if(iDw == d) {indexEigenUp[iUp] = iQMC; iUp++;continue;}
        if( eigenValUp[iUp] >= eigenValDw[iDw] ){
		    indexEigenUp[iUp] = iQMC; 
		    iUp++;
		}
		else{
		    indexEigenDw[iDw] = iQMC; 
		    iDw++;
		}
	    }
	}

    }

}
\end{hide}
\end{code}
